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Abstract. This paper is concerned with a numerical simulation of shape optimization 
in a two-dimensional viscous incompressible flow governed by Navier-Stokes equations 
with mixed boundary conditions containing the pressure. The minimization problem of 
total dissipated energy was established in the fluid domain. We derive the structures of 
shape gradient of the cost functional by using the differentiability of a minimax formu- 
lation involving a Lagrange functional with a function space parametrization technique. 
Finally a gradient type algorithm is effectively used for our problem. 
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1 Introduction 

The problem of finding the optimal design of a system for the viscous incompress- 
ible flow arises in many design problems in aerospace, automotive, hydraulic, ocean, 
structural, and wind engineering. In practise, engineers are interested in reducing the 
drag force in the wing of a plane or vehicle or in reducing the dissipation in channels, 
hydraulic values, etc. 

The optimal shape design of a body subjected to the minimum dissipate viscous 
energy has been a challenging task for a long time, and it has been investigated by 
several authors. For instance, O.Pironneau in [16, 14] computes the derivative of the 
cost functional using normal variation approach; F.Murat and J.Simon in [15] use the 
formal calculus to deduce an expression for the derivative; J.A.Bello et.al. in [1, 2, 3] 
considered this problem theoretically in the case of Navier-Stokes flow by the formal 
calculus. 
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The previous references concern Dirichlet boundary conditions for the velocity. 
However, the velocity-pressure type boundary conditions must be introduced and 
it seems more realistic in many industrial applications, such as shape optimization 
of Aorto-Coronaric bypass anastomoses in biomedical engineering. Recently, H.Yagi 
and M. Kawahara in [20] study the optimal shape design for Navier-Stokes flow with 
boundary conditions containing the pressure using a discretize-then-differentiate ap- 
proach. However, its proposed algorithm converges slowly. E.Katamine et.al. [13] 
use the differentiate-then-discretize approach with the formulae of material deriva- 
tive to study such problem involving velocity-pressure type boundary conditions with 
Reynolds number up to 100. 

In the present paper, we will use the so-called function space parametrization tech- 
nique which was advocated by M.C.Delfour and J.-P.Zolesio to solving poisson equation 
with Dirichlet and Nuemann condition (see [6]). In our paper [8, 9, 10], we apply them 
to solve a Robin problem and shape reconstruction problems for Stokes and Navier- 
Stokes flow with Dirichlet boundary condition only involving the velocity, respectively. 

However, in this paper we extend them to study the energy minimization problem 
for Navier-Stokes flow with velocity-pressure boundary conditions in despite of its lack 
of rigorous mathematical justification in case where the Lagrange formulation is not 
convex. We shall show how this theorem allows, at least formally to bypass the study 
of material derivative and obtain the expression of shape gradient for the given cost 
functional. Finally we will introduce an efficient numerical algorithm for the solution 
of such minimization problems and the numerical examples show that our proposed 
algorithm converges very fast. 

This paper is organized as follows. In section 2, we briefly recall the velocity method 
which is used for the characterization of the deformation of the shape of the domain 
and give the description of the shape minimization problem for the Navier-Stokes flow 
with mixed type boundary conditions. 

Section 3 is devoted to the computation of the shape gradient of the Lagrangian 
functional due to a minimax principle concerning the differentiability of the minimax 
formulation by function space parametrization technique. 

Finally in section 4, we give its finite element approximation and propose a gradient 
type algorithm with some numerical examples to show that our theory could be very 
useful for the practical purpose and the proposed algorithm is efficient. 

2 Preliminaries and statement of the problem 

2.1 Elements of the velocity method 

To our little knowledge, there are about three types of techniques to perform the 
domain deformation: J.Hadamard [12] 's normal variation method, the perturbation 
of the identity method by J.Simon [18] and the velocity method (see J.Cea[4] and J.- 
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P.Zolesio[6, 21]). We will use the velocity method which contains the others. In that 
purpose, we choose an open set D in R^ with the boundary dD piecewise C k , and a 
velocity space V € E k := C([0,e]; [V k (D)] N ), where e is a small positive real number 
and [V k (D)] N denotes the space of all A;— times continuous differentiable functions with 
compact support contained in D. The velocity field 

V(t)(x) = V(t,x), xeD, t>0 

belongs to [V k (D)] N for each t. It can generate transformations Tt(V)X = x(t,X) 
through the following dynamical system 

(AX 

— (t,X) = V(t,x(t)), x(0,X)=X 

with the initial value X given. We denote the "transformed domain" Tt(V)(Q) by 
Cl t (V) at t > 0, and also set its boundary T t := T t (T). 

There exists an interval I = [0, S), < S < e, and a one-to-one map Tt from D onto 
D such that 

(i) T = I; 

(ii) (t,x) ^ T t (x) belongs to C 1 (I;C k {D)) with T t {dD) = dD; 

(iii) (t,x) ^ T t ~ 1 (x) belongs to C{I;C k {D)). 

Such transformation are well studied in [6]. 

Furthermore, for sufficiently small t > 0, the Jacobian Jt is strictly positive: 

J t {x) := det|DT t (x)| = detDT t (x) > 0, 

where DTt (x) denotes the Jacobian matrix of the transformation Tt evaluated at a point 
x 6 D associated with the velocity field V. We will also use the following notation: 
DT t _1 (x) is the inverse of the matrix DTt(x) , *DT t _1 (s) is the transpose of the matrix 
DT f _1 (x). These quantities also satisfy the following lemma. 

Lemma 2.1 For any V € E k , DT t and J t are invertible. Moreover, DTt, DT t are 
in C 1 ([0,£];[C k - 1 (D)] NxN ), and J t , Jf 1 ore tnC 1 ([0,e];C*- 1 (5)). 

Now let J(O) be a real valued functional associated with any regular domain $7, we 
say that the functional J (CI) has a Eulerian derivative at £1 in the direction V if the 
limit 



Jim - [J(fl t ) - J(n)] := d J(0; V) 



exists. 

? k 



Furthermore, if the map V i— > d J(Q; V) : E — > R is linear and continuous, we say 
that J is shape differentiable at fl. In the distributional sense we have 

dJ(0;F) = (VJ,V)^ V k^) N yxV k (D) N - 
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When J has a Eulerian derivative, we say that VJ is the shape gradient of J at 0. 

Before closing this subsection, we introduce the following functional spaces which 
will be used in this paper: 



H(n) 
v g (n) 

V {Q) 



N 



{u £ H\n) 

{u G H 2 {Q) N 
{u G H 2 {Q) N 



divu = in Q, it = on T w UT S , u = g on T u }, 
u = on T w U T s , u = g on T u }, 

« = 0on r„,ur u ur s }, 



<5(0) : = J p 6 fl-i(n) : J pdx = 0{ if meas(r d ) 



0) 



2.2 Formulation of the flow optimization problem 

Let be a region of R 2 and we denote by T the boundary of 0. We suppose that is 
filled with a Newtonian incompressible viscous fluid of the kinematic viscosity v. The 
flow of such a fluid is modeled by the following system of Navier-Stokes equations, 

-diva(y,p) +By y = in 0, (2.2) 
divy = in Q, (2.3) 

where i/ denotes the velocity field, p the pressure, and a(y,p) the stress tensor defined 
by a(y,p) := —pl + 2ue(y) with the rate of deformation tensor e(y) := (Dy + *Dy)/2, 
where *Dy denotes the transpose of the matrix Dy and I denotes the identity tensor. 

Equations (2.2) and (2.3) have to be completed by the following typical boundary 
conditions: 

y = g onT u (2.4) 
y = on T s U T w (2.5) 
a(y,p) n = h on T d (2.6) 

where n denotes the unit vector of outward normal on T = T u U U T w L)T S , T u is the 
inflow boundary, the outflow boundary, r„, the boundary corresponding to the fluid 
wall and T s is the boundary which is to be optimized. We also recall that the Reynolds 
number Re is classically defined by Re = UL/u with U a characteristic velocity and L 
a characteristic length. 

For the existence and uniqueness of the solution of the Navier-Stokes system (2.2)- 
(2.6), we have the following results (see [19]). 

Theorem 2.1 Suppose that Q is of class C . For the data 

g£H 3 / 2 (D) N , [ g-nds = 0; heH^ 2 {D) N , 

there exists at least one y € H(Q) and a distribution p € L 2 (0) on Q such that (2.2)- 
(2.6) hold. Moreover, if is is sufficiently large or g and h sufficiently small, there exists 
a unique solution (y,p) G H(£l) x L 2 (ft) to the problem (2.2) -(2.6). In addition, if O 
is of class C 2 , we have (y,p) G V g (Q) x Q(Q). 
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Our goal is to optimize the shape of the domain SI which minimizes a given cost func- 
tional depending on the fluid domain and the state. The cost functional may represent a 
given objective related to specific characteristic features of the flow (e.g., the deviation 
with respect to a given target pressure, the drag, the vorticity, • • • ). 

We are interested in solving the total dissipation energy minimization problem 

min J (SI) = 2u / \e(y)\ 2 dx, (2.7) 
neo J n 

where the boundary T u U Td U T w is fixed and an example of the admissible set O is: 
O := \ SI C 1^ : T u U T d U T w is fixed, dx = constant I . 



Corollary 2.1 ([17]) Let SI be of piecewise C 1 , the minimization problem (2.7) has at 
least one solution with given area in two dimensions. 

3 Function space parametrization 

In this section we derive the structure of the shape gradient for the cost functional 
J(Sl) by function space parametrization techniques in order to bypass the usual study 
of material derivative. 

Let SI be of class C 2 , the weak formulation of (2.2)-(2.6) in mixed form is: 

' seek (y,p) G V g (Sl) x Q(Sl) such that 

< J Q [2ve(y) : e{v) +Dy y ■ v -pdivv] dx = J T h ■ vds, Vu G V (Q), (3.1) 
k f n divyqdx = 0, Vg G Q(n). 

Where in the weak form (3.1), we have used the following lemma. 

Lemma 3.1 

2 / e(y) : e(v) dx = — (Ay + V div y) ■ v dx + 2 / e(y) ■ n ■ vds. 
Jn Jn Jan 

Now we introduce the following Lagrange functional associated with (3.1) and (2.7): 

G(S},y,p,v,q) = J(Sl) - L(Sl,y,p,v,q), (3.2) 

where 

L(Sl,y,p,v,q) = / [2i/e(y) : e(v) + Dy • y ■ v — pdiv v] dx — / h vds— divyqdx. 
Jn Jv d Jn 

The minimization problem (2.7) can be put in the following form 

min min max G(Sl,y,p,v,q), (3.3) 

n&o (y,p)ev g (n)xQ(n) (v, q )<=v (n)xQ(n) 

We can use the minimax framework to avoid the study of the state derivative with 
respect to the shape of the domain. The Karusch-Kuhn- Tucker conditions will furnish 
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the shape gradient of the cost functional J(O) by using the adjoint system. Now let's 
establish the first optimality condition for the problem 

min max G(£l,y,p,v,q). (3-4) 

(y ! p)eV g (Q)xQ(Q) (U,g)eVo(0)xQ(0) 

Formally the adjoint equations are defined from the Euler-Lagrange equations of the 
Lagrange functional G. Clearly, the variation of G with respect to (v,q) can recover 
the state system (3.1). On the other hand, in order to find the adjoint state system, 
we differentiate G with respect to p in the direction 5p, 

dG f 

— (n,y,p,v,q)-5p = J Spdivvdx = 0, 

Taking 5p with compact support in Q gives 

divt> = 0. (3.5) 

Then we differentiate G with respect to y in the direction by and employ Green formula, 
dG f 

—-(£l,y,p,v, q) ■ 5y = I (-2vAy + vAv -Vq- *Dy -v + Dv-y)- bydx 
dy Jn 

— / a(v, q) ■ n ■ by ds + 4u / e(y) ■ n ■ by ds — / (y ■ n)(v ■ by) ds. 
Jdn JdVl JdfL 

Taking by with compact support in f2 gives 

- uAv + Vq + *By ■ v - Dv ■ y = -2vAy. (3.6) 
Then varying by on gives 

a(v, q) ■ n + (y ■ n)v — Ai/e(y) ■ n = 0, on T^. (3.7) 
Finally we obtain the following adjoint state system 



— div a(v, q) + *Dy ■ v — Dv ■ y = —2vAy in Q 

div v = in Q 

a(v, q) ■ n + (y ■ n)v — 4ue(y) ■ n = 0, on 

v = on T u U r„, U T s 



(3.8) 



and its variational form 

seek (v,q) <G V (tt) x Q(Q) such that V((p,ip) € V (O) x Q(fi), 

j n [2ve(v) : e(cp) + Dtp ■ y ■ v + Dy ■ <p ■ v — q div cp] dx = Au J Q e(y) : e(cp) dx, 

j n diw-^dx = 0. 

(3.9) 

We employ the velocity method to modelize the domain deformations. We only perturb 
the boundary T s and consider the mapping Tt(V), the flow of the velocity field 

V G T4d := {V € C°(0, r; C 2 (R N ) N ) : V = in the neighorhood of T U UT W U T d }. 
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We denote the perturbed domain £l t = T t (V)(Q). 

Our objective in this section is to study the derivative of j(t) with respect to t, 
where 

j(t) := min max G(flt,yt,pt,v t ,qt), (3.10) 

(y t ,Pt)ev g (n t )xQ(Q t ) (Vt,qt)ev (n t )xQ(n t ) 

(y t ,Pt) and (vt,qt) satisfy (3.1) and (3.9) on the perturbed domain Q t , respectively. 

Unfortunately, the Sobolev space V g (Clt), Vb(fit), and Q(£lt) depend on the param- 
eter t, so we need to introduce the so-called function space parametrization technique 
which consists in transporting the different quantities (such as, a cost functional) de- 
fined on the variable domain Of back into the reference domain Q which does not 
depend on the perturbation parameter t. Thus we can use differential calculus since 
the functionals involved are defined in a fixed domain with respect to the parameter 
t. 

To do this, we define the following parametrizations 

V g (n t ) = {yoT t - 1 : yeV g (n)}; 
V (n t ) = {voT t - 1 : v€V (Q)}; 
Q(n t ) = {poT t - 1 : p€Q(Q)}. 

where " o" denotes the composition of the two maps. 

Note that since T and Tf 1 are diffeomorphisms, these parametrizations can not 
change the value of the saddle point. We can rewrite (3.10) as 

j(t)= min max G(U t , yoT^^oTf 1 , voTf 1 , qoTr 1 ). (3.11) 

(y, P )ev g {n)xQ(n) (v, g )ev (n)xQ(n) 

where the Lagrangian 

G(n t , yoT t -\po Tf l ,v o T t ~\q o Tf 1 ) = h(t) + I 2 (t) + I 3 (t) 

with 

h(t) :=2u [ HyoT^dx, 

hit) := - f [2ueiv o Tf l ) : e{y o Tf 1 ) + D(y o Tf 1 ) • (y o Tf l ) • (v o Tf 1 ) 
J Sit 

- ip o Tf 1 ) div (v o Tf 1 ) - (y o Tf 1 ) ■ V(g o Tf 1 )] dx, 

and 

I 3 (t) := h-vds. 

Now we introduce the theorem concerning on the differentiability of a saddle point (or 
a minimax). To begin with, some notations are given as follows. 
Define a functional 

Q : [0, r] x X x Y -» R 
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with r > 0, and X, Y are the two topological spaces. 

For any t G [0, r], define g(t) = inf xe x sup yeY G(t, x i V) an d the sets 

X(t) = {x t G X : g(t) = sup„ ey ^(*, a;*, y)} 
^(*, x) = {y* G y : x, y*) = sup^y G(t, x, y)} 

Similarly, we can define the dual functional h(t) = sup^gy mf^x Q(t, x, y) and the 
corresponding sets 

Y(t) = {y t €Y: h(t) = inf x . eX Q(t, x, y*)} 
X(t, y) = {x t eX: Q{t, x*,y) = mi xeX g(t, x, y)} 

Furthermore, we introduce the set of saddle points 

S(t) = {(x,y) £1x7: g(t) = G(t,x,y) = h(t)} 
Now we can introduce the following theorem (see [5] or page 427 of [6]): 
Theorem 3.1 Assume that the following hypothesis hold: 
(HI) S(t)^0, t G [0, t); 

(H2) The partial derivative dtG(t,x,y) exists in [0, r] for all 



(x,y) g 



U X(t) x ^(0) 
te[o,r] 



u 



X(0) x |J 
te[o,T] 



(#5,) There exists a topology Tx on X such that for any sequence {t n : t n G [0,r]} with 
lim t n = 0, £/iere exists x° G X(0) and a subsequence {t nk }, and for each k > 1, 

n /"oo 

£/iere exists x nk G X(f nfc ) smc/i £/iat 

(%) lim x nfc = x° m i/ie 7x topology, 

n y<x> 

(it) lim inf d t G(t,x nh ,y) > d t G(0,x°,y), Vy € Y(0); 



t\,o 

(H4) There exists a topology Ty on Y such that for any sequence \t n : t n G [0, r]} mtt 
lim £ n = 0, i/iere exists y° G 7(0) and a subsequence {t nk }, and for each k > 1, 

n /"oo 

t/iere exists y nfc G Y(i nfc ) snc/i t/iai 
fi) lim y nfc = y° in i/ie TV topology, 

n /"oo 

(ii) lim sup d t G(t,x,y nk ) < d t G{0,x,y°), Vx G X(0). 

k^oc 

Then there exists (x°,y°) G X(0) x Y(0) such that 

g(t)-g(o) 



do(0) = lim : 



inf sup d t G(0, x,y) = d t G{0,x° ,y°) = sup inf d t G(0, x,y) (3.12) 



This means that {x , y ) G -X"(0) x 7(0) is a saddle point of dtG(0,x,y). 
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Following Theorem 3.1, we need to differentiate the perturbed Lagrange functional 
G(n t ,y o Tf\p o Tf\ voT t -\qo Tf 1 ). 

To perform the differentiation, we introduce the following Hadamard formula[12] 



df 



f F(t,x)dx = [ ^(t,x)dx + [ F(t, x) V ■ n t ds, (3.13) 

J fit JCit ®t JdVLt 



for a sufficiently smooth functional F : [0, r] x M. N — > R. 
By Hadamard formula (3.13), we get 

9 t G(O t ,y o r ( _1 ,po T,- 1 , v o Tf x ,q o T^ 1 ) = 4(0) + 4(0) + 4(0) + 4(0), 

where 

4(0) = Av I e(y) : e(-Dy • V ) dx + 2v [ \e(y)\ 2 V n ds; (3.14) 
Jn Jr s 

4(0) = - / [2^(-Dy • V) • e(v) + 2ue(y) ■ e(-Bv ■ V) + By ■ y ■ (-By ■ V) 
Jn 

+ B(—By -V)-yv + By (-By ■ V) ■ v - pdiv (-Bv ■ V) 
— div (—By ■ V)q — div y(—Vq ■ V) — (—Vp ■ V) div v] dx 

+ / (— 2ve(y) : e(v) — By ■ y ■ v + pdiv v + div yq)V n ds; (3.15) 

and 4(0) = 0. 

To simplify (3.14) and (3.15), we introduce the following lemma. 

Lemma 3.2 If two vector functions y and v vanish on the boundary T s and divy = 
div v = in £1, the following identities 

ByV-n = (Byn-n)V n =divyV n ; (3.16) 
e(y) : e(v) = e(y) : (e(v) ■ (n ® n)) = (e(y) ■ n) ■ (e(v) ■ n); (3.17) 
(e(y) • n) • (Bv ■ V) = (e(y) ■ n) ■ (Bv ■ n)V n = (e(y) ■ n) ■ (e(v) ■ n)V n (3.18) 

hold on the boundary T s , where the tensor product n®n:= £^fy=i n^nj. 
Using Lemma 3.1, for (3.14) we have 

4(0) = -2v ! Ay(-DyV)da? + 4i/ / (e(y ) ■ n) ■ (-By ■ V) ds + 2u I \e(y)\ 2 V n ds. 
Jn Jr s Jr s 

By the identities (3.17) and (3.18), we further get 

4(0) = -2i/ / Ay(-ByV)dx-2v [ \e(y)\ 2 V n ds. (3.19) 
Jn Jr s 
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Employing Lemma 3.1 and y\r 3 = V^jr^uruurd = 0, (3.15) can be rewritten as 

4(0) = / {(vAy -By y - Vp) ■ (-Bv • V) + divy(-Vg • V)] dx 
Jn 

+ / [(uAv + Bv ■ y - *By ■ v - Vq) ■ (-By ■ V) + d\vv(-Vp ■ V)] dx 
Jn 

- / [a(y,p) ■ n ■ (-Bv ■ V) + a(v, q) ■ n ■ (-By ■ V)} ds 
Jr 3 

— / [2ve(y) : e(v) + By ■ y ■ v — pdlvv — d\\ yq] V n ds. (3.20) 
Jys 

Since (y,p) and (v,q) satisfy (2.2)(2.3) and (3.5)(3.6) respectively, (3.20) reduces to 

4(0) = 2v f Ay(-ByV)dx 
Jn 

- [ [a(y,p) ■ n ■ (-Bv ■ V) + a(v, q) ■ n ■ (-By ■ V) + 2ve(y) : e(v)V n ] ds. (3.21) 
Jr s 

On the boundary r s , we can deduce that 

—a(y,p) ■ n ■ (—Bv ■ V) — a(v, q) ■ n ■ (—By ■ V) 

= 2u[e(y) ■ n ■ (Bv ■ V) + e(v) ■ n ■ (By ■ V)] (by (3.16)) 

= Au(e(y) ■ n) ■ (e(v) ■ n)V n (by (3.18)) 

= Aue(y):e(v)V n . (by (3.17)) 

Therefore, (3.21) becomes 

4(0) = 2v \ Ay (-By ■ V)dx + 2v / e(y) : e(v)V n ds. (3.22) 
Jn Jr s 

Adding (3.19) and (3.22) together, we finally obtain the boundary expression for the 
Eulerian derivative of J(O), 

dJ(Q; V) = 2v [ [e(y) : e(v) - \e(y)\ 2 ] V n ds, (3.23) 
Jr s 

Since the mapping V i— > d,J(£l; V) is linear and continuous, we get the expression for 
the shape gradient 

VJ = 2v[e(y) : e(v) - \e(y)\ 2 ]n (3.24) 

by (2.1). 

4 Finite element approximations and numerical Simula- 
tion 

4.1 Discretization of the optimization problem 

We suppose that is a bounded polygonal domain of ]R 2 and only consider the con- 
forming finite element approximations. Let C i? 1 (J7) Ar and Sh C L 2 (Q) be two 
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families of finite dimensional subspaces parameterized by h which tends to zero. We 
also define 

V g h ■= {u h G X h : u h = on T w U T s , u h = fir on T u }, 
V oh := {u h G X h : u h = on T w U T u U 

:= Iph € 5ft : y ph. drc = ( if meas(r d ) = 0)1 . 

Besides, the following assumptions are supposed to hold. 
(HA1) There exists C > such that for < m < I, 

mf |K-«||i < Ch m \\v\\ m+1 , Vv e H m+1 (n) N nv g (n) ■ 

Vh&v gh 

(HA2) There exists C > such that for < m < I', 

inf \\q h - q\\ < Ch m \\q\\ m , Vg G # m (ft) D Q(fi); 

(HA3) The Ladyzhenskaya-Brezzi-Babuska inf-sup condition is verified, i.e., there exists 
C > 0, such that 



inf 



J^div^dx ^ 



V h = V gh or Vnh. 



Q¥=<ih£Qho^v h eV h \\ v h\\i\\Qh\\o 



oh ■ 



The Galerkin finite element approximations of the state system (3.1) and adjoint state 
system (3.9) in mixed form are as follows 



seek {y h ,Ph) G Vgh x Qh such that V(v h ,q h ) G V oh x 
f a [2ue(y h ) : e(v h ) + T)y h ■ y h ■ v h - p h div v h ] dx = f Fd h-v h ds, 
fn divy h q h dx = 0, 



(4.1) 



and 



seek (v h ,q h ) G V oh x Q h such that V((p h ,ir h ) G V oh x 
J Q [2ue(v h ) : e(<p h ) + Dcp h ■ y h ■ v h + By h -<p h -v h - q h divcp h ]dx 

k J n divv h -K h dx = 0. 
We also have the discrete cost functional 

J h (n) = 2 [ \e(y h )\ 2 dx, 
Jn 

and the discrete shape gradient 

VJ h = 2u[e(y h ) : e(v h ) - \e(y h )\ 2 ]n 
Finally for completeness, we state the following theorem (see [11])- 



(4.2) 



(4.3) 



(4.4) 
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Theorem 4.1 Assume that the hypotheses (HAl), (HA2) and (HAS) hold. Let 

{(A, (y(X), Xp(X))); A = 1/u G A, A is a connected subsect o/M + } 

6e a branch of nonsingular solutions of the state system (3.1). Then there exists a 
neighborhood O of the origin in V g (Q) x Q(Q) and for h < ho sufficiently small a 
unique C°° branch {(A, (y h (X), Xph(X))); X G A} of nonsingular solutions of problem 
(4.1) such that 

lim sup{||y h (A) - y(A)|| 2 + \\ Ph (X) -p(A)||i} = 0. 

In addition, for the adjoint state system (3.9) and its discrete form (4.2), we have the 
similar convergence result. 

4.2 A gradient type algorithm 

For the minimization problem (2.7), we rather work with the unconstrained minimiza- 
tion problem 

min G(Q) = J (SI) + IV (Q), (4.5) 

where V(£l) := Jq dx and I is a positive Lagrange multiplier. The Eulerian derivative 
of G(n) is 

dG(n-V) = [ VG-Vds, 

where the shape gradient VG := [2ue(y) : e(v)—2v\e(y)\ 2 - l rl]n. Ignoring regularization, 
a descent direction is found by defining V = —h^S/G, and then we can update the shape 
$1 as fifc = (I + hk V)f2, where hy. is a descent step at fe-th iteration. 

However, in this article in order to avoid boundary oscillations (and irregular shapes) 
and due to the fact that the gradient type algorithm produces shape variations which 
have less regularity than the original parametrization, we change the scalar product 
with respect to which we compute a descent direction, for instance, H 1 ^) 2 . In this 
case, the descent direction is the unique element d G H 1 ^) 2 of the problem 

-Ad + d = in Q, 

< d = o, onr u ur d ur„,, (4.6) 

Dd n = -VG on IV 

To better understand the necessity of projection or smoother due to the loss of regu- 
larity, we give the following remark. 

Remark 4.1 We give a simple example to illustrate the loss of regularity. We suppose 
that the cost functional is a quadratic functional: J(x) = (Ax — b) 2 with x G H l (SX), 
A G and b G L 2 (JT2). The gradient VJ = 2(Ax - b)A G F _1 (fi) has less 

regularity than x. Then any variation using V J as the descent direction will have less 
regularity than x, therefore we need to project into H 1 ^). We refer the readers to 
see B.Mohammadi & O.Pironneau [14J and G.Dogan et.al.fl] for further discussion on 
regularity. 
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The resulting algorithm can be summarized as follows: 

(1) Choose an initial shape Qq, an initial step ho and a Lagrange multiplier Iq; 

(2) Compute the state system (3.1) and the adjoint system (3.9), then we can evaluate 
the descent direction d k by using (4.6) with f2 = £l k and I = Ik', 

(3) Set n k+1 = (I + h k d k ) Q k and l k+1 = (l k + l)/2 + e\V(Q k ) - V(Q)\/V(Q) with a 
small positive constant e, where I = — f r VJds/ f r ds and V(Q) is the given 
area of 

The choice of the descent step size h k is not an easy task. Too big, the algorithm is 
unstable; too small, the rate of convergence is insignificant. In order to refresh h k , we 
compare h k with h k -\. If (d k , d k ^i) H i is negative, we should reduce the step; on the 
other hand, if d k and d k _\ are very close, we increase the step. In addition, if reversed 
triangles are appeared when moving the mesh, we also need to reduce the step. 

In our algorithm, we do not choose any stopping criterion. A classical stopping 
criterion is to find that whether the shape gradients in some suitable norm is small 
enough. However, since we use the continuous shape gradients, it's hopeless for us to 
expect very small gradient norm because of numerical discretization errors. Instead, 
we fix the number of iterations. If it is too small, we can restart it with the previous 
final shape as the initial shape. 

4.3 Numerical results 

In all computations, the finite element discretization is effected using the Pibubble-Pi 
pair of finite element spaces on a triangular mesh, i.e., we choose the following velocity 
space Xh and pressure space S^: 

X h = {y h G (C°(n)) 2 : y h \ T G (P? T fyT G %} 
S h = {p h G C°{U) : p h \ T G Px.VT G %}, 

where 7^ denotes a standard finite element triangulation of £1, P k the space of the 
polynomials in two variables of degree < k and P^ T the subspace of P3 defined by 

Pit = {l '■ 1 = Qi + ^t, with q\ G Pi, A G R and 

4>t ^ P3, 4>T = on dT, 4>t{Gt) = 1 with Gt is the centroid of T}. 

Notice that a function like 4>t is usually called a bubble function. 

The mesh is performed by a Delaunay-Voronoi mesh generator (see [14]) and during 
the shape deformation, we utilize the a metric-based anisotropic mesh adaptation tech- 
nique where the metric can be computed automatically from the Hessian of a solution. 
We run the programs on a home PC with Intel Pentium 4 CPU 2.8 GHz and 1GB 
memory. 



13 



4.3.1 Test case 1: cannula shape optimization in Stokes flow 

We consider the shape optimization of a two-dimensional inflow cannula of a circulatory 
assist device in the biomedical applications. The geometry of the cannula f2 is depicted 
in the left picture of Figure 4.1. The boundary conditions for the problem are traction- 
free at the exit T^, no-slip at all curved walls T s , and a specified parabolic inlet velocity 
9(0, V) = ((y-2)(2.35-y),0) T 



v 




Figure 4.1: The analytic domain and the finite element mesh of the cannula. 

In this test case, we present results for two different Reynolds numbers 0.1 and 
0.01, defined as Re = d\y m \/v, where y m is the maximum velocity at the inlet T u 
and d = 0.35 is the diameter of the cannula. The domain is discretized using 448 
triangular elements and 279 nodes. Since the inertial term in (2.2) can be neglected 
when Re = 0.1 or 0.01, we can say that the blood flow in the cannula was governed by 
the Stokes equations approximately. 




Figure 4.2: The initial and optimal cannula shapes with Re = 0.1. 
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Figure 4.3: The initial and optimal cannula shapes with Re = 0.01. 



The distributions of the horizontal velocity for the initial and optimal shapes with 
Re = 0.1, 0.01 are shown in Figure 4.2 and Figure 4.3. It is clear that shape optimization 
has removed the sharp bend in the initial configuration of the cannula. 

The optimization process gave a 55.4975% reduction in the dissipated energy with 
Re = 0.1, and a 55.0392% reduction in the dissipated energy with Re = 0.01. 

4.3.2 Test case 2: Optimization of a solid body in the Navier Stokes flow 

As a second test case, we consider the isolated body problem. The schematic geometry 
of the fluid domain is described in Figure 4.4, corresponding to an external flow around 
a solid body S. We reduce the problem to a bounded domain D by introducing an 
artificial boundary dD := r n U U T w which has to be taken sufficiently far from S so 
that the corresponding flow is a good approximation of the unbounded external flow 
around S and O := D\S is the effective domain. In addition, the boundary T s := dS 
is to be optimized. 



► 9 














CD 











Figure 4.4: External flow around a solid body S. 

We choose D to be a rectangle (—0.5, 1.5) x (—0.5, 1.5) and S is to be determined 
in our simulations. The inflow velocity is assumed to be parabolic with a profile 
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g(— 0.5, y) = (0.2y 2 — 0.05, 0) T , while at the outflow boundary T^, we impose a traction- 
free boundary condition (h = 0). No-slip boundary condition are imposed at all the 
other boundaries. We further define the admissible set 

0:= {JlcK 2 : dD is fixed, the area V{tt) = 1.9} , 

which means that the target volume of S to be optimized is 0.1. 

We choose the initial shape of the body S to be a circle of center (0, 0) with radius 
r = 0.2. We present results for two different Reynolds numbers Re = 40, 200 defined 
by Re = 2r\y m \/v, where y m is the maximum velocity at the inflow T u . The finite 
element mesh used for the calculations at Re = 200 has been shown in Figure 4.5. 




Figure 4.5: the finite element mesh. 

Figure 5.6 and Figure 5.7 represent the distribution of the velocity y = (yi,y2) T 
and the pressure p for the initial shape and the optimal shape in the neighborhood of 
S for Re = 40, 200 respectively. 

We run many iterations in order to show the good convergence and stability prop- 
erties of our algorithm, however it is clear that it has converged in a smaller number of 
iterations (see Figure 5.8 and Figure 5.9). For the Reynolds numbers 40 and 200, the 
total dissipated energy reduced about 34.47% and 44.46% respectively. 

5 Conclusion 

The minimization problem of total dissipated energy in the two dimensional Navier- 
Stokes flow with mixed boundary conditions involving pressure has been presented. 
We derived the structure of shape gradient for the cost functional by function space 
parametrization technique without the usual study of the derivative of the state. Though 
for the time being this technique lacks from a rigorous mathematical framework for the 
Navier-Stokes equations, a gradient type algorithm is effectively used for the minimiza- 
tion problem for various Reynolds numbers. Further research is necessary on efficient 
implementations for time-dependent Navier-Stokes flow and much more real problems 
in the industry. 
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(a) y\ for initial shape. (b) yi for initial shape. (c) p for initial shape. 




(a) yi for initial shape. (b) 2/2 for initial shape. (c) p for initial shape. 




(d) yi for optimal shape. (e) j/2 for optimal shape. (f) p for optimal shape. 



Figure 5.7: Comparison of the initial shape and optimal shape for Re = 200. 
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Figure 5.8: Convergence history of the total dissipated 
energy for Re = 40. 




Iteration 



Figure 5.9: Convergence history of the total dissipated 
energy for Re = 200. 
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